
for i = 2:1:41
    try
        data1 = load(sprintf('%d_1.txt', i));
        data2 = load(sprintf('%d_2.txt', i));
        data3 = load(sprintf('%d_3.txt', i));
        count=count+1;
    catch
        continue; 
    end

    Fs=1000; %sampling rate
    Fc_low=0.5;% low cutoff freq
    Fc_high=10; % high cutoff freq
    order=3;
    Fn=Fs/2;
    [A,B,C,D] = cheby2(order, 30, [Fc_low Fc_high]/Fn,'bandpass');
    [filter_SOS, g] = ss2sos(A, B, C, D);
    data1_filt = filtfilt(filter_SOS, g, data1);
    [A,B,C,D] = cheby2(order, 30, [Fc_low Fc_high]/Fn,'bandpass');
    [filter_SOS, g] = ss2sos(A, B, C, D);
    data2_filt = filtfilt(filter_SOS, g, data2);
    [A,B,C,D] = cheby2(order, 30, [Fc_low Fc_high]/Fn,'bandpass');
    [filter_SOS, g] = ss2sos(A, B, C, D);
    data3_filt = filtfilt(filter_SOS, g, data3);

    figure;
    subplot(3,1,1);
    plot(data1_filt);
    title('Filtered Signal 1');
    subplot(3,1,2);
    plot(data2_filt);
    title('Filtered Signal 2');
    subplot(3,1,3);
    plot(data3_filt);
    title('Filtered Signal 3');
    
    % Adjust plot settings as needed
    xlabel('Sample');
    ylabel('Amplitude');
    sgtitle('Filtered Signals Comparison');

    t1=(0:length(data1)-1)/samplingRate;
    t2=(0:length(data2)-1)/samplingRate;
    t3=(0:length(data3)-1)/samplingRate;
    % if i==414
    %     figure;
    % 
    %     % Plot the first PPG signal
    %     subplot(2, 3, 1);
    %     plot(t1,data1*10^(-3));
    %     title('PPG Signal 1');
    % 
    %     % Plot the second PPG signal
    %     subplot(2, 3, 2);
    %     plot(t2,data2*10^(-3));
    %     title('PPG Signal 2');
    % 
    %     % Plot the third PPG signal
    %     subplot(2, 3, 3);
    %     plot(t3,data3*10^(-3));
    %     title('PPG Signal 3');
    % 
    %     % Plot the fourth PPG signal
    %     subplot(2, 3, 4);
    %     plot(t1,data1_filt*10^(-3));
    %     title('PPG Signal 4');
    % 
    %     % Plot the fifth PPG signal
    %     subplot(2, 3, 5);
    %     plot(t2,data2_filt*10^(-3));
    %     title('PPG Signal 5');
    % 
    %     % Plot the sixth PPG signal
    %     subplot(2, 3, 6);
    %     plot(t3,data3_filt*10^(-3));
    % %     hold on;
    % % 
    % %     % Mark the peaks
    % %     plot(t(peakIdx), peaks, 'ro', 'MarkerSize', 10);
    % % 
    % %     % Mark the troughs
    % %     plot(t(troughIdx), -troughs, 'bo', 'MarkerSize', 10);
    % %     title('PPG Signal 6');
    % % else
    % %     continue;
    % end
    % snr1 = snr(data1_filt);
    % snr2 = snr(data2_filt);
    % snr3 = snr(data3_filt);
    % 
    % % Find the PPG signal with the highest SNR
    % [maxSNR, maxIndex] = max([snr1,snr2,snr3]);
    % 
    % if maxIndex == 1
    %     data = data1_filt;
    % 
    %     time=t1;
    % elseif maxIndex == 2
    %     data = data2_filt;
    %     time=t2;
    % else
    %     data = data3_filt;
    %     time=t3;
    % end
    data=data1_filt;
     % Find the peaks in the PPG signal
    [peaks, locations] = findpeaks(data);
    % Getting first and second derivative

    % Calculate the first derivative of the data
    data_derivative = diff(data);
    
    % Calculate the second derivative of the data
    data_second_derivative = diff(data_derivative);
    
    % deriv_data=gradient(data);
    % deriv2_data=gradient(deriv_data);
    filter_size = 3;
    
    % Apply the moving average filter to the second derivative of the PPG signal
    smoothed_second_derivative = movmean(data_second_derivative, filter_size);
    
    %Create time vectors for the derivative signals
    t_derivative = t3(2:end);
    t_second_derivative = t_derivative(2:end);
    deriv2_copy=smoothed_second_derivative;
    % Create a new figure
    if i==419
         figure;
         %Plot the original PPG signal
         subplot(3, 1, 1);
         plot(t3, data * 10^(-3));
         title('Original PPG Signal');
     
    %     % Plot the first derivative of the PPG signal
         subplot(3, 1, 2);
         plot(t_derivative, data_derivative*10^(-3));
         title('First Derivative of PPG Signal');
     
    %     % Plot the second derivative of the PPG signal
         subplot(3, 1, 3);
         plot(t_second_derivative, smoothed_second_derivative*10^(-3));
         title('Second Derivative of PPG Signal'); 
     end

    % detecting SDPPG inflection points:
    inverted_ppg_signal = -data;
        
    [trough_values, trough_indices] = findpeaks(inverted_ppg_signal);
        
    trough1_index = trough_indices(1);
    trough2_index = trough_indices(2);

    time_cropped = t_derivative(trough1_index:trough2_index);   

    cropped_ppg_signal = data_derivative(trough1_index:trough2_index);
    [peaks, properties] = findpeaks(cropped_ppg_signal);
    t_axis = time_cropped(properties);
    
    if numel(t_axis) < 2
        length_deriv2=length(deriv2_copy);
        for p=1:length_deriv2               
          if ( deriv2_copy(p)< 0.0002 &&  deriv2_copy(p)>0) || (deriv2_copy(p)> -0.0002 &&  deriv2_copy(p)<0) %second derivative close to zero as it is discrete time series
              deriv2_copy(p)= 0;  
              % disp("here");
          end 
        end   
        inflections=[];
        numInf=1;
        for q=1:length_deriv2
            if deriv2_copy(q)==0
                inflections(numInf)=q;
                numInf=numInf+1;
            end
        end
    
        inflection_points = [];
        
        window_size = 10; % You can adjust this based on your data characteristics.
        
        for q = 1:numInf-1
            % Get the current inflection point index.
            inflection_idx = inflections(q);
            
            % Apply neighborhood smoothing around the current inflection point index.
            start_idx = max(1, inflection_idx - window_size);
            end_idx = min(length(deriv2_copy), inflection_idx + window_size);
            neighborhood = deriv2_copy(start_idx:end_idx);
            
            % Check if the neighborhood contains any valid inflection points.
            valid_indices = find(neighborhood == 0);
            if ~isempty(valid_indices)
                % Keep only one representative point in the local neighborhood.
                representative_inflection_idx = start_idx + valid_indices(1) - 1;
                inflection_points(representative_inflection_idx) = 1;
            end
        end
    
        inflection_point_indices = find(inflection_points);
    
    
    
        [peak_max,loc_max]=max(data);
        inf_point=0;
        for ij=1:numInf-1
            if(inflection_point_indices(ij)>locations(1))
                inf_point=inflection_point_indices(ij);
                break;
            end
        end
       
            AI=abs(data(inf_point)/peak_max);
      
            T_LASI=(0:length(data)-1)/samplingRate;
            LASI=1/(T_LASI(inf_point)-T_LASI(loc_max));
            disp("OTHER WAY:")
            fprintf("LASI: %f\n", LASI);

    else
        zero_crossing = find(cropped_ppg_signal > 0 & circshift(cropped_ppg_signal, -1) < 0);
        z_axis = time_cropped(zero_crossing);
        first_zero_to_second_max_time = (t_axis(2) - z_axis(1));
        
        disp("PEAKS:");
        disp(t_axis(2));
        disp("ZEROS:");
        disp(z_axis(1));
        disp(i)
        % Interpolate y-axis values corresponding to 't_axis'
        y_t_axis = interp1(time, data, t_axis);
        
        % Interpolate y-axis values corresponding to 'z_axis'
        y_z_axis = interp1(time, data, z_axis);
        
        AI=abs(y_t_axis(2)/y_z_axis(1));
        val = 1 / first_zero_to_second_max_time;
        LASI=val
        fprintf("LASI: %f\n", val);
    end
 

    %CT--- Crest Time (taken from the first derivative)-----------
    %correct----- look into it

    crest_time = z_axis(1);
    %disp('Crest Time:');
    %disp(crest_time);
    
    % catch
    %     disp(i);
    %     disp(inflections);
    %     continue;
    % end
     


    % DT-----------------------------------CORRECT BUT NOT AN AVG
    %{
    dataForTroughs=data;
    [troughs,locationsT]=findpeaks(-dataForTroughs,'MinPeakHeight',-Inf);
    t_Troughs=(0:length(dataForTroughs)-1)/samplingRate;
    n=numel(locationsT);
    loc_Trough=0;
   disp("I is:")
   disp(i)
    for l=1:1:n 
        if locationsT(l)>locations(1)
            if locationsT(l)< locations(2) && locationsT(l+1)>locations(2)
                loc_Trough=l;
            end
        end
       
    end
   
    
    % Check if at least one complete cycle is present
    if numel(locationsT) < 2
        disp('Insufficient cycles in the PPG signal.');
        return;
    end
    
    % Specify the cycle for which you want to calculate the diastolic time
    cycle_index = 1; % Modify this value to select the desired cycle
    
    % Check if the specified cycle index is valid
    if cycle_index < 1 || cycle_index > numel(locationsT)-1
        disp('Invalid cycle index.');
        return;
    end
    

    % Calculate the diastolic time for the specified cycle
    start_index = locations(cycle_index);
    end_index = locationsT(loc_Trough); 
    diastolic_time_cycle = (end_index - start_index) / samplingRate;
   
    
    % if i==15
    %     disp("start_index");
    %     disp(start_index);
    %     disp("end_index");
    %     disp(end_index);
    %     disp("diastolic time");
    %     disp(diastolic_time_cycle);
    %     disp("troughs");
    %     disp(troughs);
    %     fprintf('locationsT: %s\n', mat2str(locationsT));
    % end

    % Display the diastolic time for the specified cycle
    %disp('Diastolic Time for the Specified Cycle:');
    %disp(diastolic_time_cycle);

    %}
    % Delta T------------ (from first derivative of PPG signal)- correct
    
    % Find the zero crossings in the first derivative
    zero_crossings = find(data_derivative(1:end-1) .* data_derivative(2:end) < 0);
    
    % Check if at least 4 zero crossings are present
    if numel(zero_crossings) < 4
        disp('Insufficient zero crossings in the first derivative of the PPG signal.');
        return;
    end
    
    % Calculate the delta T between the 2nd and 4th zero crossing
    start_index = zero_crossings(2);
    end_index = zero_crossings(4);
    % disp(start_index)
    % disp(end_index)
    delta_T = (end_index - start_index) / samplingRate;
    
    % Display the delta T
    % disp('Delta T:');
    % disp(delta_T);
    
    % PEAKS AND TROUGHS OF 2ND DERIV:
    [peaks_deriv2,loc_peaks_deriv2]=findpeaks(smoothed_second_derivative, 'MinPeakHeight', 0.0005, 'MinPeakDistance', 200);
    [troughs_deriv2,loc_troughs_deriv2]=findpeaks(-smoothed_second_derivative,'MinPeakHeight', 0.0005, 'MinPeakDistance', 200);
    
    % if numel(peaks_deriv2)>6
    %    disp(i);
    % end

    % if numel(troughs_deriv2)>5
    %     disp("extra trough");
    %     disp(count);
    % end
    % if numel(peaks_deriv2)>7
    %     disp("extra peak");
    %     disp(i);
    %     if i==179
    %         figure;
    %         subplot(2, 3, 1);
    %         plot(t1,data1*10^(-3));
    %         title('PPG Signal 1');
    % 
    %         % Plot the second PPG signal
    %         subplot(2, 3, 2);
    %         plot(t2,data2*10^(-3));
    %         title('PPG Signal 2');
    % 
    %         % Plot the third PPG signal
    %         subplot(2, 3, 3);
    %         plot(t3,data3*10^(-3));
    %         title('PPG Signal 3');
    % 
    %         % Plot the fourth PPG signal
    %         subplot(2, 3, 4);
    %         plot(t1,data1_filt*10^(-3));
    %         title('PPG Signal 4');
    % 
    %         % Plot the fifth PPG signal
    %         subplot(2, 3, 5);
    %         plot(t2,data2_filt*10^(-3));
    %         title('PPG Signal 5');
    % 
    %         % Plot the sixth PPG signal
    %         subplot(2, 3, 6);
    %         plot(t3,data3_filt*10^(-3));
    %         figure;
    %         % Plot the original PPG signal
    %         subplot(3, 1, 1);
    %         plot(t3, data * 10^(-3));
    %         title('Original PPG Signal');
    % 
    %         % Plot the first derivative of the PPG signal
    %         subplot(3, 1, 2);
    %         plot(t_derivative, data_derivative*10^(-3));
    %         title('First Derivative of PPG Signal');
    % 
    %         % Plot the second derivative of the PPG signal
    %         subplot(3, 1, 3);
    %         plot(t_second_derivative, data_second_derivative*10^(-3));
    %         title('Second Derivative of PPG Signal'); 
    %     end
    % 
    % 
    % end

    % IH
    IH=0;
    
    %fprintf('Peaks: %s\n', mat2str(peaks));
    s1=numel(locations);%s1=number of peaks, s=1
   % for j= 1:1:s1-1
   %     IH=IH+peaks(1,j);
   % end
   % IH=IH/(s1-1);
    %IL:
    IL=0;
    [peaks2,locations2]=findpeaks(-data,'MinPeakHeight',-Inf);
    s3=numel(locations2);
    peaks2=abs(peaks2);
    
    %fprintf('Peaks2: %s\n', mat2str(peaks2));
    for k= 1:1:s3-1
        IL=IL+peaks2(1,k);
    end
    IL=IL/(s3-1);
    %PPG intensity ratio (PIR):
    PIR=IH/IL;
    % figure;
    %plot(t3,data*10^(-3));

    % MEU:
    fft_data=fft(data);
    fft_data(1)=0;
    AC_component=real(ifft(fft_data));
    [peaks3,locations3]=findpeaks(AC_component);
    s5=numel(locations3);
    MEU=0;
    for l= 1:1:s5-1
        MEU=MEU+peaks3(1,l);
    end
    MEU=MEU/(s5-1);
    %womersley number (alpha):
    alpha=IL*sqrt((1060*heart_rate(1,count))/MEU);
    %t=(0:length(data)-1)/samplingRate;
    t=(0:(1/samplingRate):(length(data)-1));
    SUT=0;
    count2=0;
    for m=1:1:2
        try
            if locations2(m)>locations(1)
                SUT=SUT+t(locations(m))-t(locations2(m));
                count2=count2+1;
            else
                continue;
            end
        catch 
            continue;
        end 
    end 
    SUT=abs(SUT/(count2));
    if SUT<0
        disp(i);
    end
    % if i==115
    %     figure;
    %     plot(t_Troughs,dataForTroughs*10^(-3));
    %     title("data excluding everything before first peak");
    %     figure;
    % 
        %Plot the first PPG signal
        % figure;
        % subplot(2, 3, 1);
        % plot(t1,data1*10^(-3));
        % title('PPG Signal 1');
        % 
        % % Plot the second PPG signal
        % subplot(2, 3, 2);
        % plot(t2,data2*10^(-3));
        % title('PPG Signal 2');
        % 
        % % Plot the third PPG signal
        % subplot(2, 3, 3);
        % plot(t3,data3*10^(-3));
        % title('PPG Signal 3');
        % 
        % % Plot the fourth PPG signal
        % subplot(2, 3, 4);
        % plot(t1,data1_filt*10^(-3));
        % title('PPG Signal 4');
        % 
        % % Plot the fifth PPG signal
        % subplot(2, 3, 5);
        % plot(t2,data2_filt*10^(-3));
        % title('PPG Signal 5');
        % 
        % % Plot the sixth PPG signal
        % subplot(2, 3, 6);
        % plot(t3,data3_filt*10^(-3));
    % 
    % else
    %     continue;
    % end

 if (LASI==-1.101321586)
         figure;
       fprintf("LASI of wrong: %f\n", LASI);
       disp("I is")
       disp(i)
    %     % Plot the first PPG signal
         subplot(2, 3, 1);
         plot(t1,data1_filt*10^(-3));
         title('PPG Signal 1');
    % 
    %     % Plot the second PPG signal
         subplot(2, 3, 2);
         plot(t2,data2_filt*10^(-3));
         title('PPG Signal 2');
    % 
    %     % Plot the third PPG signal
         subplot(2, 3, 3);
         plot(t3,data3_filt*10^(-3));
         title('PPG Signal 3');
    % 
    %     % Plot the fourth PPG signal
         subplot(2, 3, 4);
         plot(data*10^(-3));
         title('PPG Signal 4');
    % 
    %     % Plot the fifth PPG signal
         subplot(2, 3, 5);
         plot(time_cropped,cropped_ppg_signal*10^(-3));
         title('cropped');
    % 

    % %     hold on;

    
    end

end



